clc; clear; close all;

% 固定参数（论文给定或参考前面设置）
r = 3.99;           % Quantum Logistic 参数r
b = 30;             % Quantum Logistic 参数b
r1 = -10;           % 混合比例参数
b1 = -3;            % 偏移常数
a = 1;              % 3D-SIMM参数a
b2 = 2*pi;          % 3D-SIMM参数b2
c = 11.5;           % 3D-SIMM参数c

% 初始值（复数形式或实数）
x0 = 0.3;
y0 = 0.5;
z0 = 0.6;

sequence_lengths = round(logspace(4, 6, 50)); % 10^4到10^6的50个点
times = zeros(size(sequence_lengths));

for i = 1:length(sequence_lengths)
    N = sequence_lengths(i);
    
    tic;
    seq = generate_CTBCS_sequence(N, x0, y0, z0, r, b, r1, b1, a, b2, c);
    times(i) = toc;
end

% 绘图
figure;
plot(sequence_lengths, times, 'LineWidth', 1.5);
xlabel('The size of sequence');
ylabel('Time (s)');
title('The calculation time of CTBCS');
grid on;
set(gca, 'XScale', 'log'); % 横轴对数刻度

%% 生成CTBCS序列的函数
function seq = generate_CTBCS_sequence(N, x0, y0, z0, r, b, r1, b1, a, b2, c)
    x = x0; y = y0; z = z0;
    seq = zeros(1, N);
    for k = 1:N
        state = CTBCS(x, y, z, r, b, r1, b1, a, b2, c);
        x = state(1);
        y = state(2);
        z = state(3);
        seq(k) = real(x); % 取实部
    end
end

%% CTBCS系统迭代函数（你提供的）
function out = CTBCS(x, y, z, r, b, r1, b1, a, b2, c)
    F = QuantumLogistic(x, y, z, r, b);
    G = SIMM(x, y, z, a, b2, c);

    arg_x = pi * (r1 * F(1) + (1 - r1) * G(1)) - b1;
    arg_y = pi * (r1 * F(2) + (1 - r1) * G(2)) - b1;
    arg_z = pi * (r1 * F(3) + (1 - r1) * G(3)) - b1;

    x_new = cos(arg_x);
    y_new = cos(arg_y);
    z_new = cos(arg_z);

    out = [x_new; y_new; z_new];
end

%% QuantumLogistic函数（复共轭版本）
function out = QuantumLogistic(x, y, z, r, b)
    xd = conj(x);
    zd = conj(z);

    x_new = r * (x - abs(x)^2) - r * y;
    y_new = -y * exp(-2*b) + exp(-b) * r * ((2 - x - xd)*y - x*zd - xd*z);
    z_new = -z * exp(-2*b) + exp(-b) * r * (2*(1 - xd*z)*z - 2*x*y - x);

    out = [x_new; y_new; z_new];
end

%% SIMM函数（3D-SIMM模型）
function out = SIMM(x, y, z, a, b2, c)
    if x == 0, x = 1e-6; end
    if y == 0, y = 1e-6; end
    if z == 0, z = 1e-6; end

    x1 = sin(b2*z) * sin(c / x);
    y1 = sin(b2*x1) * sin(c / y);
    z1 = sin(b2*y1) * sin(c / z);

    x_new = a * x1;
    y_new = a * y1;
    z_new = a * z1;
    out = [x_new; y_new; z_new];
end

